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Abstract. In this paper we present a study of pattern formation in bidimensional systems 
with competing short-range attractive and long-range repulsive interactions. The interaction 
c/j ' parameters are chosen in such a way to analyse two different situations: the spontaneous 
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pattern formation due to the presence of strong competing interactions on different length 
scales and the pattern formation as a response to an external modulating potential when the 
system is close to its Lifshitz point. We compare different Monte Carlo techniques showing 
*^ , that Parallel Tempering technique represents a promising approach to study such systems 

and we present detailed results for the specific heat and the structural properties. We also 
present random phase approximation predictions about the spontaneous pattern formation (or 
microphase separation), as well as linear response theory predictions about the induced pattern 
formation due to the presence of an external modulating field. In particular we observe that 
^ . the response of our systems to external fields is much stronger compared to the response of a 

' Lennard-Jones fluid. 

o 

^ 1. Introduction 

d ; Spontaneous pattern formation can easily be observed in system far from thermal equilibrium 
3 due to energy and matter fluxes ([1] and references therein), but it can be observed also in 
'"O | equilibrium states as a result of the competition between interactions operating on different 
q | length scales [1, 2, 3, 4, 5, 6, 7]. Patterns appear both in two and three dimensional systems 
^ ; even if, in the last decades, bidimensional or quasi bidimensional systems have grown in 
. ^ | importance also due to the possibility of technological applications. In particular the possibility 
X ' of manipulating the pattern formation throughout the application of external fields has gained 
d ; great interest. Typical experiments involve fluid films subject to interfering multiple laser beams 
or ferrofluids in magnetic fields, so that the external agent can manipulate the morphology of 
clusters as well as the nature of transitions (e.g. the laser induced freezing) [8, 9, 10]. 
The competing terms are typically due to a short-range attraction, which tends to induce a 
macroscopic phase separation, and a longer-range repulsion, present in addition to the hard 
core contribution, which causes a frustration of the system favouring the formation of smaller 
clusters of particles. The short-range attraction is often explained in terms of van der Waals 
forces, while the long-range repulsion is often related to effective or actual dipolar interactions as 
in Langmuir monolayers and in magnetic films [3, 11, 12]. Other source of long-range repulsion 
is a screened Coulomb interaction of charged macroparticles ([5, 13] and references therein). 
In this paper we present continuum simulation results of a bidimensional system subject to 
an effective potential which includes competing interactions. Specific features arise if the 



Bidimensional fluid system with competing interactions 



2 



repulsion is truly long-range like the dipolar one. We do not explore this aspect since we 
are interested in generic aspects of competition, so both attraction and repulsion are assumed 
to have a finite range. At first, the potential parameters have been set to study the spontaneous 
breaking of symmetry leading to pattern formation (or microphase separation) without external 
field. Simulations have been done adopting different Monte Carlo techniques; in particular, the 
Parallel Tempering scheme has been explored. Finally the parameters of interactions have 
been chosen so to study the pattern formation induced by an external modulating potential, 
when the fluid is near its Lifshitz point. In fact it is known from previous studies [14] that 
specific features appear in the properties of the system even when the long-range repulsion is 
not strong enough to give rise to microphase separation and a normal liquid-vapor transition 
is present. Close to this Lifshitz point the coexistence curve is extremely flat and the region of 
large compressibility increases enormously compared to the normal case. Under these condition 
we might expect a very large response to an external modulating field. Both cases have been 
analysed in different thermodynamic states. 

The paper is organized as follows. In section 2 we discuss and compare different Monte 
Carlo techniques to simulate systems with competing interactions. Computational results are 
described in section 3; in section 4 we compare simulation data with theoretical predictions 
based on the random phase approximation and describe a very simple model to explain why 
we observe different patterns as density increases. Section 5 is dedicated to the analysis of the 
pattern induced by an external modulating potential. A comparison with a standard Lennard- 
Jones fluid response is also shown. In section 6 we draw conclusions. 



2. Microphases: model system and computational details 

The model system studied in this paper is similar to that in [2]. Unlike the authors [2] we 
have focused our attention on the thermodynamic and structural features of the system. We 
have analysed such quantities on a large range of temperatures, obtaining information on the 
microphase region and also on the transition among the ordered states and the homogeneous 
one. The particles interact via a spherically symmetric pairwise-additive potential which is a 
sum of two parts: 

U(r) = U sr (r) + U lr (r). (1) 
The short-range contribution U sr is the hard disk potential, while the long-range one is: 

MO = -^ ex P(-f ) + ^ ex P(-f )> ( 2 ) 

JIq -ILqi It .j, -Thy 

in which r is the interparticle distance, a is the hard disk diameter; the subscript letters a and 
r refer to the "attraction" and "repulsion", so R a , R r and e a , e r are respectively the ranges and 
the strengths of the long-range interactions. In the first part of the paper we study the case 
e a = e r , so that J °° Ui r (r) dr = 0. We have performed few computations of control for the 
ranges R a = 2a and R r = 4a as in [2] but all the results reported here have the interaction 
ranges equal to R a = la and R r = 2a, in order to reduce the computational cost of simulations 
without affecting the basic physics of the system. In the figures and in all that follows we will 
use the reduced variables T = kT/U c , U = U/U c (with U c = \Ui r (r = a)\ also referred to as 
contact value of the potential), p =< p > a 2 (with < p >= N/A, N number of particles and A 
simulation box area). 
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Figure 1. Interaction potential for a pair of particles, 
defined into the text. 



The ranges of the interactions are 



N=400 <U ex /N>=-0.1 331 (5) 
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Figure 2. Snapshots relative to the state p = 0.4 T = 0.5, using different number of particles 

in MC simulations. There are no relevant modifications to the shape of the striped pattern. 

The wave vector connected to periodicity of the pattern is the same in both cases: k p a = 0.60 



We have used different strategies of simulation. At first we performed computer simulations 
using the standard Metropolis Monte Carlo algorithm (in the following MC) [15] at a fixed 
number of particles N, area A and temperature T, with periodic boundary conditions and 
minimum image convention. We applied a spherical cutoff R c (always smaller than half of the 
simulation box side), that is we set the pair potential U(r) to zero for r > R c and then we 
included tail corrections into the computation of potential energy to compensate for the missing 
long-range part of U\ r due to the cutoff. We used N = 400 but in preliminary runs we adopted 
iV = 1600 to test the effect of the finite size of simulation box and in figure 2, for instance, 
we have plotted the snapshots relative to a state exhibiting a striped pattern. We can see that 
the shape of the pattern does not change: stripes remain parallel to each other; the excess 
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internal energy per particle is the same (within the statistical uncertainty); the wave vector 
which identify the periodicity of the pattern (that we will define more precisely in section 3) 
is the same in both cases and it corresponds to k p a ~ 0.60; changes in the radial distribution 
function computed with different number of particles are less than 2%. 

Since we are dealing with systems subject rather to long-range interactions, we have tested 
different cutoff values inside the simulations. In particular, using a cutoff too short, we can 
completely miss the features of the state, simply because the long range repulsion becomes too 
weak or it is even absent. In particular at R c = 2.5a (so that the particles do not feel at all the 
repulsive hump of Ui r ) the system displays a standard liquid- vapor phase transition when T is 
low; the same happens at R c = 3a; at R c = 5a we can observe the formation of domains whose 
shape is mostly irregular, while at R c = 10a and R c = 15a striped patterns have formed the 
shape of which is very regular (see figure 3). About the last two cases we note that an increase 
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Figure 3. Snapshots relative to the state p = 0.4 T = 0.5 with different cutoff R c . 



of the 50% of the cutoff has produced an increase of ~ 6% in the pattern period. Such feature, 
however, does not affect the overall picture that can be traced about such phases by the present 
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study, which is mainly aimed to characterize general behavior relative to the pattern formation. 
Ergodicity problems were checked by starting test runs from completely different initial 
conditions: we have used both homogeneous configurations (with particles set on a square 
lattice covering the entire simulation box) and totally segregated configurations (in which all 
the particles belong to the same liquid-like cluster). In all the cases simulations converged to 
the same state. 

The standard MC method follows a Markov chain through the configuration space. Each MC 
move consists of a single particle displacement and the move is accepted with probability 1 if the 
resulting energy change AE < 0, and accepted with probability exp(— AE/kT) otherwise. This 
method produces the correct equilibrium ensemble probability in the limit of infinite simulation 
time. Inside this scheme the passage of particles from a cluster to an other is not very frequent 
when T is very small, so that the system remains "frozen" into a certain configuration for a 
long time (i.e. for many MC steps per particle). We supplemented this conventional method 
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Figure 4. An example of a droplets pattern 



with additional MC moves designed to accelerate the motion through the configuration space. 
Thus we introduced a collective move (or clusters' move) that tries to update simultaneously 
the position of many particles. In the following we will refer to this technique as MC C ;. At this 
point we want to emphasize that, here and in the following, the term cluster will be used to 
identify a generic collection of particles irrespective of its morphology. When we intend to refer 
to a particular cluster morphology we will explicitly speak of circular domains (or droplets) 
and stripes. 

Into the MC c z scheme we choose to identify the clusters using a geometrical criterion: particles 
% and j belong to the same cluster if their separation d = |r» — ij| < R . Typically R Q = 1.5a, 
corresponding to an attractive interaction between particles. In other words, the particles 
belonging to the same cluster are connected by a percolative path. Having identified the 
clusters, we must ensure that any cluster move will not prevent finding the inverse move at a 
later time, which means that cluster moves must not change the particle connectivity. This 
is achieved by displacing the center of mass of each cluster uniformly in a box of fixed side, 
but forbidding cluster attachment to preserve detailed balance [16, 17]. Inside this scheme, 
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Figure 5. (a): instantaneous potential energy per particle; (b): number of clusters; the 

attempted cluster moves are a 20% of the total MC moves; the attempted swap moves in 
PT are ~ 0.05% of the total MC moves. MCS stands for MC steps per particle. The initial 
configuration was a homogeneous square lattice. 



in which we use the detailed balance principle and assume a symmetric transition matrix [18] 
for the probability of performing a trial move, the acceptance criterion reduces once again to 
mm(l, exp(— AE)), where AE is the energy change between the new and the old configurations. 
The collective moves represent the 20% of the total. 

In this way large changes of the system can be obtained in a single MC step. Such moves are 
favourable mainly at low temperature where the mean displacement of a cluster can be 1 — 2 
order of magnitude bigger than the single particle one. At high temperature this technique is 
not so effective because the greater probability of clusters' overlapping. This scheme enhances 
the sampling of the configuration space with respect to the simple MC technique. 
In addition, to be sure to sample sufficiently the configuration space and that we were not 
trapped into a local minimum of the free energy landscape, we adopted the Parallel Tempering 
algorithm (PT) ([18, 19] [20] and references therein). The PT technique was developed for 
dealing with the slow dynamics of disordered spin systems. The PT algorithm simultaneously 
simulates a set of M identical non interacting replicas of a system, each of them at a different 
temperature. Periodically a swap between configurations belonging to different replicas is 
attempted. The exchange is accepted in a Metropolis fashion applied to an extended ensemble 
which is a combination of subsystems that are NVTi ensembles [18]. The acceptance criterion 
is: 



w„ 



(3) 



1 A m , n < 

exp(-A„ v „) A mj „ > 

where A mtn = (/3 n — fl m )(E m — E n ), with E m E n energies of the replicas of index m n and (3 is the 
inverse of temperature. The success of PT relies on the fact that the temperature range covers 
values high enough to pass every free energy barrier. The outcome of such calculations is, in 
principle, equilibrium configurations in the canonical ensemble at the M different temperatures. 
In our simulation we took M = 31 and we covered the range 0.32 < T < 0.92. Each replica 
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evolved according to a MC c i scheme with 0.05% attempted swap moves on the total MC 
moves. Exchange events were examined only between subsystems that are nearest neighbours in 
temperature: % and i + 1 [21]. Comparing PT with the techniques previously described we find 
that it is very efficient to improve convergence and sample equilibrium configurations, specially 
for the low temperature states. See, for example, figure 5 where we plotted the number of 
clusters and the potential energy per particle for the thermodynamic state p = 0.1 and T = 0.32 
(a snapshot of which is shown in figure 4). In figure 6(a) we plotted the covariance relative to 
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Figure 6. (a): covariance versus distance in MCS (MC steps per particle) between 

configurations generated with different MC schemes; u = (U/N). Filled symbols: 1/e decay 
(read text) with respect to the initial value, (b): estimator of variance versus the length of 
blocks used into data blocking scheme for computation of statistical errors.; LL is the number 
of configurations into each block. 



U/N obtained with different simulation schemes versus the distance ( expressed in MC steps) 
between two different configurations. On each line the filled symbols mark the point where the 
covariance drops by a factor 1/e with respect to the initial value: the configurations generated 
with the PT algorithm become uncorrelated faster than they do with the other methods. In 
figure 6(b) we plotted the estimators for the variance, using the data blocking technique, for 
the computation of statistical averages and errors: if the simulation is sufficiently long the 
estimator tends toward a constant [18]. We see that for a simulation run of fixed length, the 
PT curve is the flattest one; i.e., using PT technique we need shorter runs. A typical length of 
our PT runs is ~ 5 • 10 6 MC steps per particle. 



3. Simulation results 



We report our results for the interaction potential of equation 1; the parameters of the potential 
are e a — e r — 1, while the ranges are R a = la and R r = 2a. Initial investigations were 
done using the standard MC algorithm, from which we located the region of phase diagram 
corresponding to different pattern formation. We explored the density range p = (0.1 4- 0.8) 
and the temperature range T = (0.3 -j- 1.0). At T < 0.7 we could observe the formation 
of microphases. In particular for p < 0.35 circular domains of particles (i.e. droplets) were 
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observed while for 0.35 < p < 0.6 liquid-like striped patterns formed. This is qualitatively 
in agreement with many experimental results in which, as the density increases, the cluster 
morphology changes to form pattern with lower curvature [2, 3]. At higher densities gas bubbles 
inside a liquid medium could be observed as a counterpart of the liquid droplets. The aim of the 
present computations was to characterize only the major features of the phase diagram, so that 
we cannot exclude the presence of some additional morphology mainly in the transition region 
between droplets and stripes. The characteristic dimension of the pattern is simply connected 
to the range of the potential: the droplet diameter as well as the stripe width is ~ (4-h5)<7 that 
roughly corresponds to the location of the repulsive hump of Ui r . Extensive computations were 
done with the PT and the cluster moves technique which essentially confirmed the preliminary 
runs but allowed more precise computation of the thermodynamic and structural quantities of 
the system. In this paper we report the results relative to two particular densities p = 0.1 and 
p = 0.4 in which the system exhibits droplet and striped pattern respectively. The simulation 
results shown in this section are obtained adopting the PT if not otherwise specified. 
In figure 7 we plotted the average excess internal energy per particle < U/N > for the droplet 
and stripe cases: at higher density the excess internal energy is less than it is in the low density 
case because particles, which belong to stripes at p = 0.4, feel a stronger reciprocal repulsion 
than they do when they are set into droplets at p — 0.1. In the following section we will discuss, 
within a simple model, the change of cluster morphology as the density increases. 
In order to identify the phase transition leading to the loss of order of the microphase region, 
we computed the dimensionless excess heat capacity per particle as shown in figure 8. We have 
computed in two ways C% x /Nk B at constant volume: 
CT = 1 d<U> 

Nk B Nk B dT 1 1 
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Figure 8. Upper panel: comparison of the dimcnsionless excess heat capacity per particle for 
densities p = 0.1 and p = 0.4 computed according to the equation 5; middle and bottom panels: 
comparison between the two methods of computing C v x / '{ksT) according to the equations 5 
and 4. Dashed lines are only guide to eyes. 



Since we study the system at discrete temperatures, the derivative is approximated with finite 
differences so that: 

C e v x (T n ) 1 < U(T n+l ) > - < UjT^) > 
Nk B Nk B T n+1 - r n _! 1 j 

with T n+ i > T n _i. In equilibrium these two ways of calculating C^ x should agree. In the top 
panel of figure 8 we have shown the results relative to the second method, while a comparison 
between the two methods for each density is shown in the middle and bottom panels. The 
agreement between the two methods is quite good, suggesting that the systems have well 
equilibrated at all temperatures. The peaks of the specific heat identify the phase transition 
from the microseparated region to the homogeneous one. The corresponding transition 
temperatures are T c d = 0.5 for the droplets and T c s = 0.62 for the stripes. The latter case 
shows a particularly pronounced peak and also the snapshots of the system confirm an abrupt 
passage from an ordered configuration to a disordered one (see figure 9). It is interesting to note 
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Figure 9. Snapshots relative to the microphases-liomogeneous fluid transition. 
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that the specific heat for p = 0.4 has a maximum at T = 0.34 and this seems to be connected to 
the freezing transition inside the stripes, as we will discuss below treating structural quantities 
such as the radial distribution function g(r) and the static structure factor S(k). 
In figure 10 we have plotted g(r) (computed averaging over all the directions and normalized to 
the mean density of the system) for different temperatures; in the same figure we have indicated 
the coordination number n c : 

n c = 2np / rg(r)dr, (7) 
Jo 

in which R m in is the position of the first minimum. The curves are shifted for clarity; we can 
identify two regimes: a short-range modulation due to the local structure of the fluid and a 
longer-range modulation due to the microphases formation. The latter manifests itself with a 
shallow minimum around p ~ 6cr at low temperature, which is doomed to disappear as the 
system becomes homogeneous at larger T. As temperature decreases, the second peak of g(r) 
develops a shoulder which eventually grows till to a really peak splitting. The appearance 
of this shoulder has been referred to as a structural freezing precursor both in two and three 
dimensions [22]. The shoulder appears at T ~ (0.38-7-0.40) for droplets and at T ~ (0.44 -=-0.46) 
for stripes, while the splitting is clearly visible at T = 0.32 for droplets and T = 0.38 for stripes. 
The peak splitting is consistent with a triangular lattice of period a ~ 1.05 a. 
The static structure factor has been calculated by explicit evaluation of the expression: 

N N 

S(k) = N- 1 < (J2 cos(k • r,)) 2 + sin(k ■ r,)) 2 >, (8) 

i i 

with k along different directions so, in the following, we will refer to S x (k) and S y (k) as the 
static structure factor with k along the x and y directions respectively. Due to the finite size of 
the simulation box, we can compute such quantity only for wave vectors multiple of the smallest 
vector k Q = 2n/Lb (Lb box simulation side). In figure 11 we have reported the results at one 
temperature for the two studied densities. We can recognize two regimes. For ka < 3 the peaks 
are strictly connected to the modulation of the pattern which the clusters arrange themselves 
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Figure 10. Left: g(r) for p = 0.4; right: g(r) for p = 0.1; next to the curves the temperature 
and (into parenthesis) the coordination number are indicated; the curves are shifted for clarity. 
The critical temperature is also reported. 
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Figure 11. Static structure factor at T = 0.36 computed along x and y directions for different 

densities. The insets show the structure factor compared with an estimator of the Fourier 

transform of the density profile along the same direction. 



on. We refer to the peaks at low k as to the modulation peaks. The height of the modulation 
peak strongly depends on the direction of the stripes and the direction of k along which we 
are computing S(k), i.e., if the stripes are close to one of the principle axis of the simulation 
box or not. For instance, at p = 0.4 plotted in figure 11, the x direction lays perpendicular 
to the stripes, while the y direction is quite parallel to them. Into the inset of the same figure 
we have also shown a comparison between S(k) and the square of the Fourier transform of 
the density modulation | < p k > 2 \/N, to underline that the peaks to a large extent are due 
just to the density modulation. Of course on a very long run the pattern will fluctuate giving 
a uniform density so < p k > will vanish whereas S(k) maintains the modulation peaks. The 
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location of the main modulation peak corresponds to a modulation period A = 2n/k ~ Ilex. 
For ka > 3, instead, the behavior of S(k) is determined mainly by the structure of the fluid 
inside the clusters. As a matter of fact, at low temperature we can observe a flattening of the 
peaks at ka > 10 till to a true peak splitting (at the lowest temperatures), which is just the 
counterpart into the Fourier space of the peaks splitting observed into the radial distribution 
function. 

Even if the interparticle interaction potential is characterized by a spherical symmetry, we 
see that, in the microphase region of the phase diagram, the pattern shape can show a breaking 
of symmetry. Such feature is particularly evident in the stripe case where the particle domains 
are aligned along a fixed direction, even if such direction is selected at random. To make clear 
this aspect we have shown in figure 12 the structure factor computed on a grid of k vectors 
(the mesh of which depends on the simulation box sides as previously explained), using the 
MC C /, at different temperatures. Below the critical temperature (upper panels) the pattern 
shape has a two-fold symmetry as we can argue from the presence of a modulation peak (and 
its specular one) only along a single direction. At T = 0.32 secondary peaks can also be seen 
due to the harmonics. The height of the modulation peaks decreases with temperature, their 
width, instead, is limited by the finite size of the simulation cell. At the temperatures higher 
than the critical one (bottom panels), the structure factor shape changes: modulation peaks 
appear in different direction till to form a ring at small k. The presence of the ring means that, 
also at high temperatures, enhanced density fluctuations are present and their characteristic 
length scale is linked to the inverse of the radius of such ring according to the formula A ~ 27r/fc. 
Some anisotropy exhibited by the height of the structure factor along the ring is still present at 
T >T C and this can be interpreted as a measure of the persistence of such density fluctuations, 
persistence which can be very long (data shown in figure 12 are, in fact, obtained via simulation 
whose run length is 5 10 6 MC steps). Similar results are obtained with PT technique. 



4. Random phase approximation and simple model for phase behavior 

The static structure factor S(k) gives many important details on the structure of a fluid; S(k) 
defined as: 

S(h) = l + ph(h), (9) 

where h(r) is the total correlation function, can be also expressed, via the Ornstein-Zernicke 
equation as 

s( *> = T=W (10) 

in which the direct correlation function c(k) appears. Equation 10 is more useful because, for 
potentials with steep repulsion plus a weak tail, we can use the following decomposition: 

c(k) = c sr (k) + ci r (k), (11) 

making a distinction between a short-range contribution (c sr ) due to the hard disk potential 
and a long-range one due to Uy. For c sr we adopted the expression given in [23], which comes 
from simulation data fitting. For the long-range contribution we have used the random phase 
approximation (RPA) q r ~ —Ui T (k)/kBT. In figure 13 we compare this approximate structure 
factor for few thermodynamic states at high temperature with simulation data. The theoretical 
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T=0.32 T=0.60 




Figure 12. Structure factor at different temperatures for p = 0.4 adopting the MC C ; technique. 
Upper panels shows temperatures below the critical one, while the bottom panels are at 
temperatures higher than T c s . 

as well as the simulation S(k) exhibit a large peak for ko < 1 which is connected at high 
temperature to a characteristic length scale of fluctuations which at low temperature becomes 
the characteristic wave vector of the patterns. The insets show a comparison between Chd and Q r : 
for ko < 1 ci r has a much stronger dependence on k than cm, so that the structure is mainly 
dominated by the long-range potential (see also [24]). Simulations however show stronger 
density fluctuations compared with the RPA results. The agreement between simulation and 
RPA about the modulation peak is better at higher temperature, since the approximation of a 
weak potential tail is better fulfilled. Within the RPA, the locus of points where S(k) diverges 
is the spinodal that we reported in figure 14, where we also reported the critical temperature 
estimated for p = 0.1 and p = 0.4 from the specific heat. The spinodal curve is quite similar to 
that of a classical fluid which undergoes only a fluid-fluid separation and it exhibits a critical 
density around p ~ 0.3. From our preliminary simulation we found a transition from circular 
domains to stripe domains just around p ~ (0.3 -j- 0.35), but we have not yet explored in 
detail this transition region. At low T due to competition between long-range repulsion and 
short-range attraction a fluid-fluid separation is forbidden, but particles arrange themselves 
on smaller clusters to decrease reciprocal repulsion. Increasing density, however, the cluster- 
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Figure 13. RPA predictions about the structure factor compared with simulation data for 

different densities and temperatures in the homogeneous regime; in the insets we show a 

comparison between the hard disk direct correlation function and the contribution due to 

long-range term of the potential. 
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Figure 14. Spinodal curve obtained studying the divergence of the structure factor according 
to equations 10 and 11. Crosses are simulation results obtained from the study of specific heat 
(see previous section). 



cluster repulsive interaction is not negligible anymore and above a certain density a new pattern 
becomes favoured in order to minimize repulsion, leading to a morphology change that in our 
case is a passage from droplets to stripes. Let us assume that the transition is dominated by 
energetic effects and we introduce a very simplified model in order to predict the stability of the 
two patterns. We treat separately the interactions among particles which belong to the same 
cluster and the interactions among particles belonging to different clusters, in order to enlighten 
the role of reciprocal interaction among clusters. The first kind of interaction is referred to with 
the subscript letter c, while the latter ones with the subscript letters cc. The excess internal 
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energy per particle can be so expressed as: 

U ex /N = Uc + x-u m (12) 

where x = Ndz/2 (N c i number of clusters and z coordination number of the clusters, i.e. 6 
for droplets supposed on a triangular lattice and 2 for stripes supposed on a regularly spaced 
grid) is proportional to the number of first neighbours clusters. Our goal is to compute such 
quantity in the density range p(0.1 -j- 0.5) for each pattern (droplets and stripes). 
Once chosen the pattern we suppose the clusters are identical, that is all the droplets have 
the same diameter and all the stripes have the same width. We also suppose that the average 
distance d among clusters is fixed and set equal to the value obtained by RPA. Since N (total 
number of particles), A simulation box area and d average distance among clusters are fixed, 
we can obtain, for each pattern, the number of clusters N c i, the number of particle n and the 
density p c i inside each cluster and their area S c [. 
An approximate expression for u c and u cc are: 



u c = (U ex /n) c = ^ Pd I U(r)g(r)dr (13) 



u. 



{U ex /n) cc = ^l [ t/MsiV^drMr 2 (14) 

&cl JS\ JS 2 , 



where U(r) is the interaction potential adopted of equation 1, the integral is over the cluster 
domain and the superscript numbers 1 e 2 underline the fact that we refer to two different 
clusters. The expressions for u c and u cc are obtained in analogy of the excess internal energy 
that can be computed for a homogeneous fluid subject to a pair-wise potential. To be very 
simple we neglect correlations, i.e. we put 

0(r)~l , , 

(2)/ i 2\ 1 y LtJ ) 

Remembering that U ex = — 9X q® (Q partition function) and that the free energy is 
A = — In Q//3, we can obtain an estimation of the free energy as: 

since the excess energy does not depend on temperature in our approximation. f} is a reference 
state with respect to which we compute the change of free energy. Assuming that the reference 
state is the ideal one, such approximation leads to A ex = U ex , so that what we will say for the 
excess internal energy can be immediately related to the free energy. 

If in equation 12 we neglect the contribution u cc due to interaction among clusters we find 
that the circular domain should be preferred for all densities (see figure 15); on the contrary 
including u cc we see that the results are totally different: there exist different density regions in 
which droplet patterns and striped patterns appear respectively preferred (figure 15), i.e. there 
exist a "critical density p D " such that for p < p the system manages to minimize its excess 
internal energy (that, in our model, is equal to the free energy) setting particles on a droplet 
pattern, while for p > p Q it achieves such energy minimization setting particles on a striped 
pattern. Inside this very simplified model the passage between droplets and stripes happens 
around p Q ~ 0.33, consistent with our simulation evidence. 
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Figure 15. Model predictions about the droplets-stripes transition, with and without cluster- 
cluster repulsive interaction. 



5. Microphases induced by external modulating potentials 

If the intensity e r of the repulsive long-range potential in equation 2 is small the system displays 
a standard liquid-vapor phase transition with a critical temperature which is a decreasing 
function of e r . The value of e r at which a fluid subject to competing interactions stops to undergo 
a standard liquid-vapor transition favouring a liquid-modulated phase is called Lifshitz point 
(LP) [14]. In RPA this point corresponds to a change of convexity of the Fourier transform of 
the interaction potential at zero wave vector. Theoretical studies of 3D systems [25] emphasize 
that, approaching LP (but without entering the microphase region yet), the inverse of the 
compressibility of the fluid has a very small value over an extended region in density and 
in temperature. That means that the system manages to achieve large density fluctuations 
without undergoing a phase transitions. Such a behavior suggests that just in this region the 
presence of an external modulation should greatly affect the system. In this section we discuss 
such situation in bidimensional systems. The interaction potential is the same as in equation 
1: the ranges R a and R r have not been changed, while the intensity of the short-range and 
long-range interactions are e a = 0.679 and e r = 0.22. Such parameters correspond to the LP. 
In figure 16 we show the interaction potential and its Fourier transform for different values of 
A = e r /e a . We added, along the x direction, an external modulating potential of the form : 



in which the strength e m is in units of the contact value U c (defined in section 2), X{ refers 
to the i-particle abscissa; the period of the function is P = A/2 that, in the present case, we 
have taken equal to half of the simulation box side. Simulations have been done adopting MC C | 
scheme (described in section 2) in a canonical ensemble and with N = 400 particles. 
Adopting the random phase approximation described in section 4 we estimate the critical point 
position at p c ~ 0.28 and T c ~ 0.95. From simulation we have verified that the fluid, without 
modulating potential, undergoes a standard liquid- vapor transition below the critical point. 
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Figure 16. Interaction potential with different value of the parameters relative to the strength 

of the attractive and repulsive component (left) and its Fourier transform (right). The ranges 
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Figure 17. Radial distribution function for our system at the thermodynamic state p = 0.28 
T = 1.1 next LP compared with a typical behavior of the same quantity for a classical LJ fluid 
(without external modulation). 



We studied a temperature range rather close to the estimated T c along the critical isochore. 
Without external modulation and T > T c the fluid phase is characterized by a radial distribution 
function very slowly decaying to unity, as we can see from figure 17 where we plotted g(r) for 
the thermodynamic state p = 0.28 T = 1.1 in comparison with a case relative to a bidimensional 
Lennard- Jones system (LJ) confirming the great influence of the density fluctuations in this 
particular region. Switching on the modulation, particles tend to organize in stripes located in 
correspondence of the minimum of the external potential (see figure 18 for a snapshot of the 
system and a typical density profile along the modulation direction). Obviously a stripe pattern 
becomes better defined as e m increases. In figure 19 an example of structure factor is reported 
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Figure 18. Left: snapshot relative to our system for p = 0.28 T = 1.1 and e m = 0.1. Right: 
number or particles versus position along the modulation axis normalized to the mean number 
of particles without external potential. 
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Figure 19. Static structure factor along x and y directions without external modulation (left 

panel) and with external modulation (right panel); the position of the modulation peak is 

indicated. Dotted lines are only a guide to eyes. 



with and without external modulation. For e m = x and y directions are equivalent since the 
fluid is homogeneous and isotropic; the large values of S(k) at small k are due to the proximity 
of the critical point. For e m = 0.1 at ka ~ 0.33 (corresponding to the modulation wave vector of 
the external potential) S(k) has particularly increased with respect to the counterpart computed 
along the y axis. 

To test the effect of the density fluctuations close to the LP, we have compared our system with 
the bidimensional LJ, studying the response to an external modulation at different temperatures 
and for various intensities of the external potential. To simplify the comparisons we have also 
expressed temperatures and modulation strength in terms of the critical temperature for each 
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Figure 20. Upper panels: LRT predictions and simulation results about the Fourier transform 
of the density profile along the modulation axis in correspondence at the modulation wave vector 
for a LJ fluid (left) and for our system (right), versus reduced temperature t = (T — T c )/T c at 
fixed e* n . Bottom panel: LRT and simulation results versus e* m at fixed temperature. For each 
system temperatures and strengths of the modulating potential are expressed in units of their 
own critical temperature. 



system. LJ data are collected in table 1, where T = k B T/w and e m = e m /w (w depth well), 
t = (T — T C LJ )/T C LJ , e* m = e m w/T^ J ; estimation of the critical point for a LJ system in 2D can 
be found in literature, according to which p^ J ~ 0.35 and T c iJ ~ 0.47 [26] (for a critical review 
about 2D LJ see also [27, 28]. Our system data are collected in table 2 where t = (T — T c )/T c 
and e* m = e m U c /T c . 

From simulations we computed the Fourier transform of the density profile along the modulation 
direction (< pk > s ' m ). Then we applied the linear response theory (LRT) according to which 

< p k >=< p k > -f3pS(k) U m {k), (18) 
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Table 1. Temperature and external modulation strength, used in simulations for the LJ system, 
in units of the well depth and in units of the LJ critical temperature. See the text for explicit 
definitions. 



T 


t 




e* 

m 


0.507 


0.079 


0.05 


0.106 


0.519 


0.104 


0.10 


0.213 


0.546 


0.162 


0.15 


0.319 


0.596 


0.268 






0.800 


0.702 







Table 2. Temperature and external modulation strength, used in simulations using our 
competing interaction potential, in units of both the contact value U c and the critical 
temperature. See the text for explicit definitions. 



T 


t 




e* 


1.020 


0.074 


0.1 


0.023 


1.050 


0.105 


0.2 


0.046 


1.100 


0.158 


0.4 


0.091 


1.200 


0.263 







where the subscript o denotes a canonical average over the unperturbed system and U m (k) is 
the Fourier component of the external potential at wave vector k; since the latter corresponds to 
a homogeneous fluid phase we can set < pk > = 0. The static structure factor which appears in 
the equation 18 has been estimated via simulation at e m = (having verified that < pk > s im = 
within the errors). In the upper panels of figure 20 we plotted the results obtained for LJ and 
for our system relative to the modulation wave vector as a function of temperature. In both 
cases we see that, for small modulation strength e m , LRT is a very good approximation even 
close to the critical point. In the bottom panel of figure 20 we plotted a comparison between 
LRT predictions and simulation data as a function of e* m for a given t: the response to the 
external potential of the system with competing interaction is stronger than it is in the LJ 
system and the deviation from the linear theory starts at a smaller value of e* m . For example 
in our case at e* m = 0.091 < pk > pea k — 83 and LRT is already not very accurate, while in 
the LJ system at e* m = 0.106 < pk > pea k — 20 and LRT still works very well. Obviously as 
increases too much LRT is not more accurate also for the LJ system. 

6. Conclusions 

In this work we have presented detailed simulation data relative to the study of spontaneous 
pattern formation in bidimensional system with competing short-range attractive and long- 
range repulsive, showing that the PT technique is a promising approach, since this technique 
reaches equilibrium configurations faster than other algorithms and it improves the statistics 
quality. At the densities of our computations, we have observed the formation of droplet-like 
or stripe-like patterns. The location of the transition temperature from the modulated region 
to the homogeneous phase has been clearly identified through the position of the peak in the 
specific heat. 

We have also analysed how the arising of patterns affects structural quantities such as g(r) and 
S(k). The radial distribution function exhibits a long-range modulation due to the underlying 
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pattern, which is doomed to disappear as temperature increases. The main feature, into the 
static structure factor, is the appearance of peaks (modulation peaks) at small wave vector. In 
particular for p = 0.4 at T < T c S(k), computed on the grid of small wave vectors, is anisotropic 
in k and has a two-fold symmetry due to the breaking of symmetry that the system undergoes 
to, while at T > T c S(k) is isotropic with a unique strong peak at small wave vector. This 
signals that strong density fluctuations are present at temperature above T c where no patterns 
are present. Besides, in the striped case we observe in the S(k) profile, at low temperatures, 
the presence of many peaks at wave vectors which are multiple of the modulation period. The 
origin of such harmonics can be traced back to the rather abrupt density variation between its 
maximum and minimum values. This effect could be checked through scattering experiments. 
Some harmonics peaks are also present at p = 0.1 in the droplet phase but these peaks are 
weaker and wider. 

In a small system, like ours, the direction of the striped pattern appears to be stabilized by 
the periodic boundary conditions; on the contrary in a very large system we expect that the 
direction of the stripes will fluctuate in time. It should be an important test determining the 
time scale of such fluctuations since they might influence properties like the birefrigence that 
could be experimentally measured. 

For the homogeneous state the random phase approximation works quite well in giving the 
position of the fluctuation peak and it is semiquantitative for its amplitude. Even if RPA 
predicts only the existence of a region in which the fluid is unstable with respect to a density 
modulation, without giving any information on the pattern shape, it is interesting to note the 
presence of a critical point in the spinodal curve given by RPA and the critical density is into 
the region in which, via simulation, we observed the passage from droplets to stripes. At last 
a very simple model, to mimic the droplet-stripe transition, has been described. The most 
relevant aspect is that the shape of the pattern is an expedient to optimize and reduce cluster 
reciprocal repulsion. The development of a microscopic theory able to describe the modulated 
phases is an important open issue. 

Finally we have studied the fluid behavior when it is close to its Lifshitz point but it is still in a 
homogeneous phase which it is strongly affected by density fluctuations. We have so analysed 
the system response to the presence of an external agent, such as a modulating potential. 
Comparing the results with those deriving from a standard LJ fluid, subject to the same 
external modulation, we have seen that our system is much more affected by the action of the 
field. In other words, next to the Lifshitz point the fluid experiments large density fluctuations 
which can be easily driven into a stable microphase by the presence of a very weak external 
potential. 
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